!*******************************************************************************
!** Main code to write the read_data information for lammps run with alkenes
!** If more than 1 molecule will be involved in the system, molecules are 
!** randomly rotated and translated to fill the space given by your box; the
!** distance between C--C, C--H, H--H from different molecules are guaranteed
!** to be greater than 1.55, 1.20, 1.0 Angstrom, respectively.
!*******************************************************************************
program main
use prec
use cell
use general
use iounits
use wallinfo
use regioninfo
use waterinfo
use H2MolInfo
implicit none
   !----------------------------------------------------------------------------
   integer             :: i, moltyp = 1, istr, iend, iprev, ijob = 1, imol, nsect
   integer             :: putasis, narg
   character (len=256) :: input
   logical             :: logdum, loverlap, skip_overlap_4_mol = .False.
   real(q)             :: rndnum(6), zlo, zhi, inbox(3)
   !----------------------------------------------------------------------------
   subname = 'main'
   info    = 'Wrong input!'
   call RANDOM_SEED( )
   !
   narg = iargc()
   i = 1
   if (i.le.narg) then
      call getarg(i, input)
      i = i + 1
      select case ( trim(input) )
      case ('-s', '-S', '--skip-overlap-check')
           skip_overlap_4_mol = .True.
      case ('-h', '-H', '--help')
           call help()
           stop
      end select
   endif
   !
   write( stdout, 10 )
   write( stdout, 20 ) "Please select the molecule you want:"
   write( stdout, 25 ) "1. Ethylene (C2H4);"
   write( stdout, 25 ) "2. Butene (C4H8);            12. Hydrogenated hexene dimer;"
   write( stdout, 25 ) "3. Hexene (C6H12);           13. Beta hexene;"
   write( stdout, 25 ) "4. Octene (C8H16);           14. Gamma hexene;"
   write( stdout, 25 ) "5. Decene (C10H20);          15. Dye molecule;"
   write( stdout, 25 ) "6. Dodecene (C12H24);        16. Mol Read from file;"
   write( stdout, 25 ) "7. LC 5CB;                   17. H2;"
   write( stdout, 25 ) "8. H2O;"
   write( stdout, 25 ) "9. Hexane (C6H14);"
   write( stdout, 25 ) "0. Exit."
   write( stdout, 26 ) "Your choice [1]: "
   read ( stdin, '(A)', iostat=ioerr ) input
   write( stdout, 10 )
   call error( subname, info, ioerr )
   !
   if ( input.ne.'' ) then
      read( input, *, iostat=ioerr ) moltyp
      call error( subname, info, ioerr )
   endif
   if ( moltyp.eq.0 ) stop
   !
   select case (moltyp)
   case ( 1:6 )
      call alkene_init( moltyp*2 )
   case ( 7 )
      call lc5cb_init()
   case ( 8 )
      call h2o_init()
   case ( 9, 12 )
      call alkane_init( (moltyp-6)*2 )
   case ( 13, 14 )
      call beta_alkene_init( moltyp-12 )
   case ( 15 )
      call dye_init();
   case ( 16 )
      call ReadMol_init();
   case ( 17 )
      call H2_Init();
   case default
      stop 'Unsupport molecular type!'
   end select
   !
   call disp_mol_info
   !
   write( stdout, 10  )
   write( stdout, 20 ) "Please select the job to do:"
   write( stdout, 25 ) "1. Generate data file;"
   write( stdout, 25 ) "2. Calculate density;"
   write( stdout, 25 ) "3. Calculate mass;"
   write( stdout, 25 ) "4. Display pair coeffs;"
   write( stdout, 25 ) "0. Exit."
   write( stdout, 26 ) "Your choice [1]:"
   read ( stdin, '(A)', iostat=ioerr ) input
   write( stdout, 10 )
   call error( subname, info, ioerr )
   !
   if ( input.ne.'' ) then
      read( input, *, iostat=ioerr ) ijob
      call error( subname, info, ioerr )
   endif
   select case ( ijob )
   case ( 1 ) ! go on

   case ( 2 ) ! calculate density
      write( stdout, 26 ) 'How many moleculars do you have ?[1]: '
      read ( stdin, '(A)', iostat=ioerr ) input
      call error( subname, info, ioerr )
      nmols = 1
      if ( input.ne.'' ) then
         read( input, *, iostat=ioerr ) nmols
         call error( subname, info, ioerr )
      endif
      if ( nmols.lt.1 ) then
         info = 'Too few moleculars!'
         call error( subname, info, 1)
      endif
      mass    = dble(nmols) * molmass
      delv    = 0.D0
      write( stdout, 26 ) 'Please input your box size or volume in Angs: '
      read ( stdin, '(A)', iostat=ioerr ) input
      call error( subname, info, ioerr )
      read ( input, *, iostat=ioerr ) box
      if ( ioerr.ne.0 ) then
         read( input, *, iostat=ioerr ) volume, delv
         if ( ioerr.ne.0 ) then
            read(input, *, iostat=ioerr ) volume
            call error( subname, info, ioerr )
         endif
      else
         volume  = box(1)*box(2)*box(3)
      endif
      density = mass/NA/volume*10.D0
      write( stdout, 60 ) density
      if ( delv.gt.0.D0 ) write( stdout, 62 ) density * delv / volume
      write( stdout, 65 ) den0
      write( stdout, 70 ) (density-den0)/den0*100.D0
      stop
   case ( 3 ) ! calculate mass
      write( stdout, 26 ) 'How many moleculars do you have ?[1]: '
      read ( stdin, '(A)', iostat=ioerr ) input
      call error( subname, info, ioerr )
      nmols = 1
      if ( input.ne.'' ) then
         read( input, *, iostat=ioerr ) nmols
         call error( subname, info, ioerr )
      endif
      if ( nmols.lt.1 ) then
         info = 'Too few moleculars!'
         call error( subname, info, 1)
      endif
      mass    = dble(nmols) * molmass
      write( stdout, 75 ) mass
      write( stdout, 76 ) mass / NA
      stop
   case ( 4 )
      write( stdout, 10 )
      do i = 1, ntypeMol
         write( stdout, 80 ) i, i, epison(i), sigma(i)
      enddo
      i = ntypeMol + 1
      write( stdout, 80 ) i, i, 0.10D0, 4.05D0
      write( stdout, 10 )
      stop
   case default ! stop
      stop
   end select
   !
   write( stdout, 26 ) 'How many moleculars do you want to have? [1]: '
   read ( stdin, '(A)', iostat=ioerr ) input
   call error( subname, info, ioerr )
   nmols = 1
   if ( input.ne.'' ) then
      read( input, *, iostat=ioerr ) nmols
      call error( subname, info, ioerr )
   endif
   if ( nmols.lt.1 ) then
      info = 'Too few moleculars required!'
      call error( subname, info, 1)
   endif
   !
   mass    = dble(nmols) * molmass
   box     = (mass/NA/den0*10.D0)**(1.D0/3.D0)
   natom   = natomMol  * nmols
   nbonds  = nbondsMol * nmols
   nangls  = nanglsMol * nmols
   ndihrs  = ndihrsMol * nmols
   !
   natmtyp = ntypeMol
   nbndtyp = nbndtypMol
   nangtyp = nangtypMol
   ndihtyp = ndihtypMol
   !
   ! Confirm if a wall will be added on top and bottom of the box.
   call init_wall()
   !
   if ( LAddWalls ) then
      nX      = int( box(1) / xLen + 0.5D0 )
      nY      = int( dble(nX) * xLen / yLen + 0.5D0 )
      !
      write( stdout, 100 ) nX, nY
      write( stdout, 110 )
      read( stdin, '(A)', iostat=ioerr ) input
      call error( subname, info, ioerr )
      if ( input.ne.'' ) then
         read(input, *, iostat=ioerr) nX, nY
         call error( subname, info, ioerr )
      endif
      if ( nX.lt.1.or.nY.lt.1 ) then
         info = 'The dimension of the wall must be greater than 0!'
         call error( subname, info, ioerr )
      endif
      !
      box(3)  = (box(1)*box(2)*box(3)) / (dble(nX*nY)*xLen*yLen)
      box(1)  = dble(nX) * xLen
      box(2)  = dble(nY) * yLen
   endif
   logdum  = any( box.gt.minlen )
   !
   ! Ask for the box size; at least one dimension should be greater
   ! than the length of the selected molecule.
   do while ( .true. )
      if ( LAddWalls ) then
         write( stdout, 120 ) dble(nX) * xLen, dble(nY) * yLen
      endif
      if ( logdum ) then
         write( stdout, 30 ) "Recommended box sizes are : ", box
         write( stdout, 35 ) "Please input your desired box sizes [", box, "]:"
      else
         write( stdout, 27 ) minlen
         write( stdout, 26 ) "Please input your desired box sizes (3F):"
      endif
      read( stdin, '(A)', iostat=ioerr ) input
      call error( subname, info, ioerr )
      if ( input.ne.'' ) then
         read(input, *, iostat=ioerr) inbox
         if (ioerr.ne.0) then
            read(input,*,iostat=ioerr) inbox(3)
            call error( subname, info, ioerr )
            inbox(1:2) = box(1:2)
         endif
         box = inbox
      endif
      if ( any(box.le.0.D0) ) then
         write(stdout, 20 ) "Box size must be greater than 0! Try again."
         cycle
      endif
      if ( any(box.gt.minlen) ) exit
   enddo
   !
   xLen    = box(1)/dble(nX)
   yLen    = box(2)/dble(nY)
   volume  = box(1)*box(2)*box(3)
   density = mass/NA/volume*10.D0
   rbox    = 1.D0 / box
   !
   ! To create the position of the walls
   zlo = 0.D0
   zhi = box(3)
   !
   if ( LAddWalls ) then
      call add_wall( zlo, zhi )
      write(stdout, 210)
      read(stdin, '(A)', iostat=ioerr) input
      if (ioerr.eq.0.and.(input.eq.'y'.or.input.eq.'Y') ) call writemap(natom)
   endif
   !
   zboxlow  = zlo
   zboxhigh = zhi
   !
   call disp_sys_info
   !
   ! Ask if molecules will be added in a limited region
   call init_region()
   !
   nsect = max(nmols/10,2)
   allocate( atpos(3, natom), attyp(natom), atchg(natom), molID(natom), stat=ioerr)
   if ( nbonds.gt.0 ) allocate( bonds(3, nbonds) )
   if ( nangls.gt.0 ) allocate( angls(4, nangls) )
   if ( ndihrs.gt.0 ) allocate( dihrs(5, ndihrs) )
   !
   info    = 'Not enough memory to allocate the necessary variables.'
   call error( subname, info, ioerr )
   !
   imol = 0
   !
   ! If only one molecule is required, ask if put as is or randomly
   if ( nmols.eq.1 ) then
      write(stdout, 28)
      read(stdin, '(A)', iostat=ioerr) input
      if ( ioerr.ne.0 ) then
         putasis = 0
      elseif ( input.eq.'' ) then
         putasis = 0
      else
         read(input,*,iostat=ioerr) putasis
         if (ioerr.ne.0) putasis = 0
      endif
   endif
   !
   write( stdout, 33 )
   !
   if ( nmols.gt. 1.or.putasis.eq.0 ) then ! put molecules randomly rotated and translated inside the box.
      do i = 1, nmols
         !
         istr     = (i-1) * natomMol + 1
         iend     = i     * natomMol
         iprev    = istr  - 1
         ! To rotate and translate the molecule randomly, make sure that no atom
         ! will overlap.
         do while ( .true. )
            loverlap = .false.
            call RANDOM_NUMBER( rndnum )
            rndnum(1:3) = rndnum(1:3) * TPI
            call rotate_molecular   ( rndnum(1), rndnum(2), rndnum(3), atposMol, rotatedMol, natomMol)
            call translate_molecular( rndnum(4), rndnum(5), rndnum(6), rotatedMol, natomMol)
            if ( lregion ) loverlap = check_out_region(rotatedMol)
            if (LAddWalls.and.(.not.loverlap)) then
               loverlap = any(rotatedMol(3,:).lt.zboxlow).or.any(rotatedMol(3,:).gt.zboxhigh)
            endif
            if ((.not.loverlap).and.(.not.skip_overlap_4_mol)) call check_distance( iprev, loverlap )
            if (.not.loverlap) exit
         enddo
         ! position, atom type and chage of the i-th molecule
         atpos(:, istr:iend)  = rotatedMol
         attyp(istr:iend)     = attypMol
         atchg(istr:iend)     = chargeMol
         molID(istr:iend)     = imol + molIDMol
         imol = imol + nMol
         ! bonds of the i-th molecule
         if ( nbondsMol.gt.0 ) then
            istr = (i-1) * nbondsMol + 1
            iend = i     * nbondsMol
            bonds(1, istr:iend)  = bondsMol(1,:)
            bonds(2:3,istr:iend) = bondsMol(2:3,:) + iprev
         endif
         ! angles of the i-th molecule
         if ( nanglsMol.gt.0 ) then
            istr = (i-1) * nanglsMol + 1
            iend = i     * nanglsMol
            angls(1, istr:iend)  = anglsMol(1,:)
            angls(2:4,istr:iend) = anglsMol(2:4,:) + iprev
         endif
         ! dihedrals of the i-th molecule
         if ( ndihrsMol.gt.0 ) then
            istr = (i-1) * ndihrsMol + 1
            iend = i     * ndihrsMol
            dihrs(1, istr:iend)  = dihrsMol(1,:)
            dihrs(2:5,istr:iend) = dihrsMol(2:5,:) + iprev
         endif
         !
         if (mod(i, nsect).eq.0) write(stdout,'(".",$)')
      enddo
   else   ! Put one molecular as is.
      atpos = atposMol
      attyp = attypMol
      atchg = chargeMol
      if ( nbondsMol.gt.0 ) bonds = bondsMol
      if ( nanglsMol.gt.0 ) angls = anglsMol
      if ( ndihrsMol.gt.0 ) dihrs = dihrsMol
      molID = molIDMol
   endif
   write( stdout, 34 )
   !
   ! Insert water molecules if necessary
   call water_init()
   if ( lwater ) call add_water()
   !
   ! Insert H2 molecules if necessary
   call H2Mol_Init()
   if ( LH2Mol ) call Add_H2Mol()
   !
   outfile = "data."//trim(molname)
   write( stdout, 26 ) "Please input the output file name ["//trim(outfile)//"]:"
   read ( stdin, '(A)', iostat=ioerr ) input
   call error( subname, info, ioerr )
   if ( input.ne.'' ) then
      read( input, *, iostat=ioerr ) outfile
      call error( subname, info, ioerr )
   endif
   ! To write the data file
   write( title, 40 ) nmols, trim(molname), trim(WallTitle)
   call output
   ! Write the final configuration in xyz format as well.
   write( title, 50 ) nmols, trim(molname), trim(WallTitle)
   call writexyz
   !
10 format(10x, 80("="))
20 format(10x, A )
25 format(15x, A )
26 format(10x, A, $ )
27 format(10x, "The maximum length of your molecule is: ", F10.4)
28 format(/,10x,"Only one molecule is required, put it randomly (0) or as is (1) ? [0]: ",$)
30 format(10x, A, F10.4," x",F10.4," x", F10.4 )
33 format(/,10x, "Now to generate the system, wait for a while", $ )
34 format(" Done!")
35 format(10x, A, 3F10.4, A, $ )
40 format("# Data file for ", I8, 1x, A, 1x, "molecules",A )
50 format("System of",I8,1x,A," molecules",1x,A )
60 format(10x, "The density of your system is   : ", F12.6, " g/cm^3" )
62 format(10x, "The variation of density   is   : ", F12.6, " g/cm^3" )
65 format(10x, "The experimental value is       : ", F12.6, " g/cm^3" )
70 format(10x, "Difference to experimental value: ", F12.6, " %"      )
75 format(10x, "Total mass of your system is    : ", G12.6, " g/mol"  )
76 format(43x, F13.6, " x 10^(-26) kg" )
80 format(10x, "pair_coeff", 1x,I2,1x,I2,2(3xF10.4) )
   !
100 format(10x, "The suitable dimension of the wall might be:",I3," x",I3 )
110 format(10x, "If this is not what you want, input yours  : ", $       )
120 format(10x, "The X-Y dimensions confined by the wall are: ", 2F10.4  )
210 format(10x, "Do you want to output the FFT map of wall atoms? (y/n)[n]:", $)
   !
stop
   !
contains
   !
   subroutine help()
   write(*, *) "  alkene"
   write(*, *) "Code to generate the data file for molecule systems"
   write(*, *) "  Usage:"
   write(*, *) "    alkene [options]"
   write(*, *) "Options available:"
   write(*, *) "  -s   Skip overlap check while inserting molecules. Default: not set."
   write(*, *) "  -h   Print this help."
   write(*, *)
   end subroutine
end program
